library(elevatr)
library(movecost)
library(sf)
library(rgdal)
library(leastcostpath)
library(sp)
library(raster)
library(tidyverse)
library(ggplot2)
library(dplyr)
library(terra)
library(ncdf4) 
library(exactextractr)  
library(rgeos)
library(PBSmapping)
library(geosphere)

setwd("Wilfahrt_APSR Replication Data")

### SPLIT GRID BY COUNTRY -> CALCULATE AREA TO ASSIGN GRID COUNTRY
{
  Country <- shapefile("Data/Grid Covariates/Admin Units/SampledCountries.shp")
  
  sf_use_s2(FALSE)
  
  Grid_10km <- shapefile("Data/Grid Shapefiles/Grid_10km.shp")
  Grid_10km$centroids <- getSpPPolygonsLabptSlots(Grid_10km)
  
  Country2 = st_as_sf(Country)
  Grid_10km2 = st_as_sf(Grid_10km)
  Atlas_8hr2 = st_as_sf(Atlas_8hr)
  
  Grid_10km2$Grid_area <- st_area(Grid_10km2) #area of each "square"
  Grid_10km2$area_km2 <- units::set_units( Grid_10km2$Grid_area, "km^2")
 
  Country2$Country_area <- st_area(Country2) #area of each "square"
  Country2$Countryarea_km2 <- units::set_units(Country2$Country_area, "km^2")
  
  Atlas_split <- st_intersection(Country2, Atlas_8hr2)
  
  Atlas_split$Polity_8hr_country_area <- st_area(Atlas_split)
  Atlas_split$Polity_8hr_country_area_km2 <-  units::set_units(Atlas_split$Polity_8hr_country_area, "km^2")
 
  Grid_split <- st_intersection(Grid_10km2, Country2)
  Grid_split$Grid_country_area <- st_area(Grid_split) #area of each "square"
  Grid_split$Grid_country_area_km2 <- units::set_units(Grid_split$Grid_country_area, "km^2")
  
  Grid_split2 <- st_join(Grid_split, Atlas_split, join = st_intersects)

  d <- as.data.frame(Grid_split2) 
  d <- d[ -c(22) ]
  
  write.table(d, file="Data/Calculated Output/Country_Grid_Split.csv",na = "", sep=",",row.names=F)
}

## SET UP GRID  -----------------------------------------------------------
{
Grid_10km <- shapefile("Data/Grid Shapefiles/Grid_10km.shp")
Grid_10km$centroids <- getSpPPolygonsLabptSlots(Grid_10km)

Centroid_num <- Grid_10km$centroids
centroid.sp<-Centroid_num
colnames(centroid.sp) <- c("longitude", "latitude")
df <- data.frame(centroid.sp)
centroid_points <- st_as_sf(df, coords = c("longitude","latitude"), crs = 4326)   

Grid_10km$Grid_area_km <- area(Grid_10km) /1000000

Grid_0_5degree <- shapefile("Data/Grid Shapefiles/Grid_0_55_degrees.shp")
Grid_1degree <- shapefile("Data/Grid Shapefiles/Grid_1_degree.shp")
Grid_2degree <- shapefile("Data/Grid Shapefiles/Grid_2_degree.shp")
Grid_3degree <- shapefile("Data/Grid Shapefiles/Grid_3_degree.shp")
Grid_5degree <- shapefile("Data/Grid Shapefiles/Grid_5_degree.shp")

Grid_1degree_2 <- intersect(Grid_1degree, Grid_10km)
Grid_1degree_2$Grid_1deg_Area_km <- area(Grid_1degree_2) /1000000
DF <- as.data.frame(Grid_1degree_2)
write.table(DF, file="Data/Calculated Output/Grid_1_degree_MERGE.csv", na = "", sep=",",row.names=F)

Grid_2degree_2 <- intersect(Grid_2degree, Grid_10km)
Grid_2degree_2$Grid_2deg_Area_km <- area(Grid_2degree_2) /1000000
write.table(DF, file="Data/Calculated Output/Grid_2_degree_MERGE.csv", na = "", sep=",",row.names=F)

Grid_05degree_2 <- intersect(Grid_0_5degree, Grid_10km)
Grid_05degree_2$Grid_0_5deg_Area_km <- area(Grid_05degree_2) /1000000
DF <- as.data.frame(Grid_05degree_2)
write.table(DF, file="Data/Calculated Output/Grid_0_5_degree_MERGE.csv", na = "", sep=",",row.names=F)
}

## IMPORT KINGDOM ESTIMATES -----------------------------------------------------------
{
#Atlas data - all polities
  Atlas_6hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Atlas_6hr.shp")
  Atlas_8hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Atlas_8hr.shp")
  Atlas_10hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Atlas_10hr.shp")
  Atlas_6hr$Atlas_6hr_area_km <- area(Atlas_6hr) /1000000
  names(Atlas_6hr)[names(Atlas_6hr) == "Polity"] <- "Polity_6hr"
  Atlas_8hr$Atlas_8hr_area_km <- area(Atlas_8hr) /1000000
  names(Atlas_8hr)[names(Atlas_8hr) == "Polity"] <- "Polity_8hr"
  Atlas_10hr$Atlas_10hr_area_km <- area(Atlas_10hr) /1000000
  names(Atlas_10hr)[names(Atlas_10hr) == "Polity"] <- "Polity_10hr"
  
#Atlas data -  composite polities (Sokoto, Adamawa, Mossi, Yoruba)
  Atlas_Comp_6hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Composite_Polities/Atlas_Composite_6hr.shp")
  Atlas_Comp_8hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Composite_Polities/Atlas_Composite_8hr.shp")
  Atlas_Comp_10hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Composite_Polities/Atlas_Composite_10hr.shp")
  Atlas_Comp_6hr$Atlas_Comp_6hr_area_km <- area(Atlas_Comp_6hr) /1000000
  names(Atlas_Comp_6hr)[names(Atlas_Comp_6hr) == "Polity"] <- "Atlas_Comp_6hr"
  Atlas_Comp_8hr$Atlas_Comp_8hr_area_km <- area(Atlas_Comp_8hr) /1000000
  names(Atlas_Comp_8hr)[names(Atlas_Comp_8hr) == "Polity"] <- "Atlas_Comp_8hr"
  Atlas_Comp_10hr$Atlas_Comp_10hr_area_km <- area(Atlas_Comp_10hr) /1000000
  names(Atlas_Comp_10hr)[names(Atlas_Comp_10hr) == "Polity"] <- "Atlas_Comp_10hr"
  
#Atlas data - Polities late 1800s
  Atlas_1880_6hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Polities_late1800s/Atlas_1880_6hr.shp")
  Atlas_1880_8hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Polities_late1800s/Atlas_1880_8hr.shp")
  Atlas_1880_10hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Polities_late1800s/Atlas_1880_10hr.shp")
  Atlas_1880_6hr$Atlas_1880_6hr_area_km <- area(Atlas_1880_6hr) /1000000
  names(Atlas_1880_6hr)[names(Atlas_1880_6hr) == "Polity"] <- "Atlas_1880_6hr"
  Atlas_1880_8hr$Atlas_1880_8hr_area_km <- area(Atlas_1880_8hr) /1000000
  names(Atlas_1880_8hr)[names(Atlas_1880_8hr) == "Polity"] <- "Atlas_1880_8hr"
  Atlas_1880_10hr$Atlas_1880_10hr_area_km <- area(Atlas_1880_10hr) /1000000
  names(Atlas_1880_10hr)[names(Atlas_1880_10hr) == "Polity"] <- "Atlas_1880_10hr"
  
#Atlas data -Polities in late 1800s with composite polities (Sokoto, Adamawa, Mossi, Maravi)
  Atlas_Comp_1880_6hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Composite_Polities_late1800s/Atlas_Composite_1880_6hr.shp")
  Atlas_Comp_1880_8hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Composite_Polities_late1800s/Atlas_Composite_1880_8hr.shp")
  Atlas_Comp_1880_10hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Composite_Polities_late1800s/Atlas_Composite_1880_10hr.shp")
  Atlas_Comp_1880_6hr$Atlas_Comp_1880_6hr_area_km <- area(Atlas_Comp_1880_6hr) /1000000
  names(Atlas_Comp_1880_6hr)[names(Atlas_Comp_1880_6hr) == "Polity"] <- "Atlas_Comp_1880_6hr"
  Atlas_Comp_1880_8hr$Atlas_Comp_1880_8hr_area_km <- area(Atlas_Comp_1880_8hr) /1000000
  names(Atlas_Comp_1880_8hr)[names(Atlas_Comp_1880_8hr) == "Polity"] <- "Atlas_Comp_1880_8hr"
  Atlas_Comp_1880_10hr$Atlas_Comp_1880_10hr_area_km <- area(Atlas_Comp_1880_10hr) /1000000
  names(Atlas_Comp_1880_10hr)[names(Atlas_Comp_1880_10hr) == "Polity"] <- "Atlas_Comp_1880_10hr"
  
# Intersected data - contested territory
  Contested_1880_8hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Contested Territory Measures/Contested_Territory_1880.shp")
  Contested_Composite_1880_8hr <- shapefile("Data/Atlas of 19th c Polities/Final Measures/Contested Territory Measures/Contested_Territory_Comp_1880.shp")
  Contested_1880_8hr$Contested_1880_8hr_area_km <- area(Contested_1880_8hr) /1000000
  names(Contested_1880_8hr)[names(Contested_1880_8hr) == "Polity"] <- "Polity_Contested_1880_8hr"
  Contested_Composite_1880_8hr$Contested_Comp_1880_8hr_area_km <- area(Contested_Composite_1880_8hr) /1000000
  names(Contested_Composite_1880_8hr)[names(Contested_Composite_1880_8hr) == "Polity"] <- "Polity_Contested_Comp_1880_8hr"
  
# Uncertainty Estimates
  UncertaintyEst_Max <- shapefile("Data/Atlas of 19th c Polities/Uncertainty Estimates/Atlas_UncertaintyEstimates_Maximum.shp")
  UncertaintyEst_Max$UncertaintyMax_area_km <- area(UncertaintyEst_Max) /1000000
  names(UncertaintyEst_Max)[names(UncertaintyEst_Max) == "Polity"] <- "Polity_UncertaintyMax"
  UncertaintyEst_Min <- shapefile("Data/Atlas of 19th c Polities/Uncertainty Estimates/Atlas_UncertaintyEstimates_Minimum.shp")
  UncertaintyEst_Min$UncertaintyMin_area_km <- area(UncertaintyEst_Min) /1000000
  names(UncertaintyEst_Min)[names(UncertaintyEst_Min) == "Polity"] <- "Polity_UncertaintyMin"
}

## IMPORT DVs  -----------------------------------------------------------
Nightlight_2013 <- raster("Data/Outcome Variables/Nightlight/F182013.v4/F182013.v4c_web.avg_vis.tif")
Nightlight_1992 <- raster("Data/Outcome Variables/Nightlight/F101992.v4/F101992.v4b_web.avg_vis.tif")
EWI <- shapefile("Data/Outcome Variables/EstimatedWealth/EstimatedWealthData.shp")



## IMPORT COVARIATES -----------------------------------------------------------
{
Diamonds <- shapefile("Data/Grid Covariates/Diamonds/DIADATA Data/DIADATA.shp")
Gems <- shapefile("Data/Grid Covariates/Gem Data/GEMDATA.shp")
Copper <- shapefile("Data/Grid Covariates/USGS Minerals/USGS_Copper.shp")
Iron <- shapefile("Data/Grid Covariates/USGS Minerals/USGS_IRon.shp")
Tin <- shapefile("Data/Grid Covariates/USGS Minerals/USGS_Tin.shp")
Gold <- shapefile("Data/Grid Covariates/USGS Minerals/USGS_Gold.shp")
Silver <- shapefile("Data/Grid Covariates/USGS Minerals/USGS_Silver.shp")
Coastline <- shapefile("Data/Grid Covariates/Water/Africa_Coastline2.shp")
Water <- shapefile("Data/Grid Covariates/Water/Africa_waterbody.shp")
MajorRivers <- shapefile("Data/Grid Covariates/Water/mrb_rivers.shp")
AfricanCapitals <- shapefile("Data/Grid Covariates/Admin Units/African_Capitals.shp")
AfricanBorders <- shapefile("Data/Grid Covariates/Admin Units/African Borders2.shp")
SlaveExports_Nunn <- shapefile("Data/Grid Covariates/Nunn/Murdock_Nunn Merge.shp")
WhitesVeg <- shapefile("Data/Grid Covariates/Whites Vegetation/WhitesVegetation.shp")

Water2 <- intersect(Water, Grid_10km)
Water2$Water_Area_km <- area(Water2) /1000000
DF <- as.data.frame(Water2)
df2 <- DF %>% group_by(Grid_ID) %>% summarize(Water_Area_km2 = sum(Water_Area_km))
require(sp)
Grid_10km <- merge(Grid_10km,df2, by="Grid_ID")
Grid_10km$Grid_PercH20 <- (Grid_10km$Water_Area_km2/Grid_10km$Grid_area_km)*100

LandSuitability <- nc_open("Data/Grid Covariates/Land Suitability/raw_data_land_suit_land_suit_0.50x0.50.nc")
print(LandSuitability)
LandSuit_Lon <- ncvar_get(LandSuitability, "longitude")
LandSuit_Lat <- ncvar_get(LandSuitability, "latitude", verbose = F)
LandSuit_t <- ncvar_get(LandSuitability, "time")
LandSuit <- ncvar_get(LandSuitability, "data")
fillvalue <- ncatt_get(LandSuitability, "data", "_FillValue")
fillvalue
LandSuit[LandSuit == fillvalue$value] <- NA
LandSuitability <- raster(t(LandSuit), xmn=min(LandSuit_Lon), xmx=max(LandSuit_Lon), ymn=min(LandSuit_Lat), ymx=max(LandSuit_Lat), crs = 4326)
LandSuitability_2 <- crop(LandSuitability, Grid_10km, snap="out")

Precipitation <- nc_open("Data/Grid Covariates/Precipitation/chirps-v2.0.annual.nc")
Precipitation_lon <- ncvar_get(Precipitation, "longitude")
Precipitation_lat <- ncvar_get(Precipitation, "latitude", verbose = F)
Prec.array <- ncvar_get(Precipitation, "precip") 
dim(Prec.array) 
fillvalue <- ncatt_get(Precipitation, "precip", "_FillValue")
fillvalue
Prec.array[Prec.array == fillvalue$value] <- NA
Prec.array <- Prec.array[, , 1] 
Precipitation <- raster(t(Prec.array), xmn=min(Precipitation_lon), xmx=max(Precipitation_lon), ymn=min(Precipitation_lat), ymx=max(Precipitation_lat), crs = 4326)
Precipitation2  <- terra::flip(Precipitation, direction="y")
Precipitation_2 <- crop(Precipitation2, Grid_10km, snap="out")
}

## SAMPLE BY GRID  -----------------------------------------------------------
# Polity estimates
{
Atlas_6hr_Intersect <- intersect(Atlas_6hr, Grid_10km)
Atlas_8hr_Intersect <- intersect(Atlas_8hr, Grid_10km)
Atlas_10hr_Intersect <- intersect(Atlas_10hr, Grid_10km)
Atlas_6hr_Intersect$Atlas_6hr_area_km_Intersect <- area(Atlas_6hr_Intersect) /1000000
Atlas_8hr_Intersect$Atlas_8hr_area_km_Intersect <- area(Atlas_8hr_Intersect) /1000000
Atlas_10hr_Intersect$Atlas_10hr_area_km_Intersect <- area(Atlas_10hr_Intersect) /1000000

Atlas_Comp_6hr_Intersect <- intersect(Atlas_Comp_6hr, Grid_10km)
Atlas_Comp_8hr_Intersect <- intersect(Atlas_Comp_8hr, Grid_10km)
Atlas_Comp_10hr_Intersect <- intersect(Atlas_Comp_10hr, Grid_10km)
Atlas_Comp_6hr_Intersect$Atlas_Comp_6hr_area_km_Intersect <- area(Atlas_Comp_6hr_Intersect) /1000000
Atlas_Comp_8hr_Intersect$Atlas_Comp_8hr_area_km_Intersect <- area(Atlas_Comp_8hr_Intersect) /1000000
Atlas_Comp_10hr_Intersect$Atlas_Comp_10hr_area_km_Intersect <- area(Atlas_Comp_10hr_Intersect) /1000000

Atlas_1880_6hr_Intersect <- intersect(Atlas_1880_6hr, Grid_10km)
Atlas_1880_8hr_Intersect <- intersect(Atlas_1880_8hr, Grid_10km)
Atlas_1880_10hr_Intersect<- intersect(Atlas_1880_10hr, Grid_10km)
Atlas_1880_6hr_Intersect$Atlas_1880_6hr_area_km_Intersect <- area(Atlas_1880_6hr_Intersect) /1000000
Atlas_1880_8hr_Intersect$Atlas_1880_8hr_area_km_Intersect <- area(Atlas_1880_8hr_Intersect) /1000000
Atlas_1880_10hr_Intersect$Atlas_1880_10hr_area_km_Intersect <- area(Atlas_1880_10hr_Intersect) /1000000

Atlas_Comp_1880_6hr_Intersect <- intersect(Atlas_Comp_1880_6hr, Grid_10km)
Atlas_Comp_1880_8hr_Intersect <- intersect(Atlas_Comp_1880_8hr, Grid_10km)
Atlas_Comp_1880_10hr_Intersect <- intersect(Atlas_Comp_1880_10hr, Grid_10km)
Atlas_Comp_1880_6hr_Intersect$Atlas_Comp_1880_6hr_area_km_Intersect <- area(Atlas_Comp_1880_6hr_Intersect) /1000000
Atlas_Comp_1880_8hr_Intersect$Atlas_Comp_1880_8hr_area_km_Intersect <- area(Atlas_Comp_1880_8hr_Intersect) /1000000
Atlas_Comp_1880_10hr_Intersect$Atlas_Comp_1880_10hr_area_km_Intersect <- area(Atlas_Comp_1880_10hr_Intersect) /1000000

Contested_1880_8hr_Intersect <- intersect(Contested_1880_8hr, Grid_10km)
Contested_Composite_1880_8hr_Intersect <- intersect(Contested_Composite_1880_8hr, Grid_10km)
Contested_1880_8hr_Intersect$Contested_1880_8hr_area_km_Intersect <- area(Contested_1880_8hr_Intersect) /1000000
Contested_Composite_1880_8hr_Intersect$Contested_Composite_1880_8hr_area_km_Intersect <- area(Contested_Composite_1880_8hr_Intersect) /1000000

UncertaintyEst_Max_Intersect <- intersect(UncertaintyEst_Max, Grid_10km)
UncertaintyEst_Min_Intersect <- intersect(UncertaintyEst_Min, Grid_10km)
UncertaintyEst_Max_Intersect$UncertaintyEst_Max_area_km_Intersect <- area(UncertaintyEst_Max_Intersect) /1000000
UncertaintyEst_Min_Intersect$UncertaintyEst_Min_area_km_Intersect <- area(UncertaintyEst_Min_Intersect) /1000000
}

# Outcome & Control Variables
{
Grid_10km$Nightlight_20132_mean <- exact_extract(Nightlight_2013,Grid_10km,fun="mean")
Grid_10km$Nightlight_1992_mean <- exact_extract(Nightlight_1992,Grid_10km,fun="mean")
EWI_Intersect <- intersect(EWI, Grid_10km)

# Covariates
Grid_10km$Diamonds <- over(SpatialPolygons(Grid_10km@polygons), SpatialPoints(Diamonds))
Grid_10km$Gems <- over(SpatialPolygons(Grid_10km@polygons), SpatialPoints(Gems))
Grid_10km$Silver <- over(SpatialPolygons(Grid_10km@polygons), SpatialPoints(Silver))
Grid_10km$Copper <- over(SpatialPolygons(Grid_10km@polygons), SpatialPoints(Copper))
Grid_10km$Tin <- over(SpatialPolygons(Grid_10km@polygons), SpatialPoints(Tin))
Grid_10km$Iron <- over(SpatialPolygons(Grid_10km@polygons), SpatialPoints(Iron))
Grid_10km$Gold <- over(SpatialPolygons(Grid_10km@polygons), SpatialPoints(Gold))

Mountains_Binary <- raster("Data/Grid Covariates/Global Mountain Explorer/GlobalMountainsK3Binary/k3binary.tif")
Mountains_Binary2 <- crop(Mountains_Binary, Grid_10km, snap="out")
Mountains_Type <- raster("Data/Grid Covariates/Global Mountain Explorer/GlobalMountainsK3Classes/k3classes.tif")
Mountains_Type2 <- crop(Mountains_Type, Grid_10km, snap="out")
Malaria_2000 <- raster("Data/Grid Covariates/Malaria Area Project/202206_Global_Pf_Parasite_Rate_2000/202206_Global_Pf_Parasite_Rate_2000.tif")
Malaria_2000_2 <- crop(Malaria_2000, Grid_10km, snap="out")

PopDensity_0AD <- raster("Data/Grid Covariates/HYDE/Hyde3_2/0AD_pop/popd_0AD.asc")
PopDensity_0AD_2 <- crop(PopDensity_0AD, Grid_10km, snap="out")

PopCount_2010AD <- raster("Data/Grid Covariates/HYDE/Hyde3_2/2010_pop/popc_2010AD.asc")
PopCount_2010AD_2 <- crop(PopCount_2010AD, Grid_10km, snap="out")

Elevation <- raster("Data/Grid Covariates/Elevation/Africa_dem_15s.tif")
Ruggedness <- raster("Data/Grid Covariates/Ruggedness/tri.tif")
Ruggedness_MeanSlope <- raster("Data/Grid Covariates/Ruggedness/slope.tif")
Ruggedness_CellArea <- raster("Data/Grid Covariates/Ruggedness/cellarea.tif")

Grid_10km$mt_binary_mean <- exact_extract(Mountains_Binary2,Grid_10km,fun="mean")
Grid_10km$Grid_MtType <- exact_extract(Mountains_Type2,Grid_10km,fun="mean")
Grid_10km$Malaria_2000 <- exact_extract(Malaria_2000_2,Grid_10km,fun="mean")
Grid_10km$PopDensity_0AD <- exact_extract(PopDensity_0AD_2,Grid_10km,fun="mean")
Grid_10km$PopCount_2010AD <- exact_extract(PopCount_2010AD_2,Grid_10km,fun="mean")

Grid_10km$Elevation <- exact_extract(Elevation,Grid_10km,fun="mean")
Grid_10km$LandSuitability <- exact_extract(LandSuitability_2,Grid_10km,fun="mean")
Grid_10km$Precipitation <- exact_extract(Precipitation_2,Grid_10km,fun="mean")

Ruggedness_SlopeWeight<- overlay(Ruggedness_MeanSlope,
                                 Ruggedness_CellArea,
                               fun=function(r1, r2){return(r1*r2)})

Ruggedness_RuggedWeight<- overlay(Ruggedness,
                                  Ruggedness_CellArea,
                                 fun=function(r1, r2){return(r1*r2)})

Grid_10km$Rugged_CellAreaTtl <- exact_extract(Ruggedness_CellArea,Grid_10km,fun="sum")

Grid_10km$Rugged_SlopeTtl <- exact_extract(Ruggedness_SlopeWeight,Grid_10km,fun="sum")
Grid_10km$Rugged_TriTtl <- exact_extract(Ruggedness_RuggedWeight,Grid_10km,fun="sum")
Grid_10km$CellArea <- exact_extract(Ruggedness_CellArea,Grid_10km,fun="sum")

## NEXT CALC GRID AREA + WATER AREA
Ruggedness2 <- crop(Ruggedness, Grid_10km, snap="out")
Ruggedness_CellArea2 <- crop(Ruggedness_CellArea, Grid_10km, snap="out")
Ruggedness_MeanSlope2 <- crop(Ruggedness_MeanSlope, Grid_10km, snap="out")

WhitesEco_Intersect <- intersect(WhitesVeg, Grid_10km)
WhitesEco_Intersect$WhiteEco_area <- area(WhitesEco_Intersect)

sf_use_s2(FALSE)

Coastline2 = st_as_sf(Coastline)
WhitesVeg_Lines2 = st_as_sf(WhitesVeg_Lines)
WhitesEco_Intersect2 = st_as_sf(WhitesEco_Intersect)

centroid_points$d_coast_m <- st_distance(x = centroid_points, y = Coastline2)
centroid_points$d_ecoborder_m <- st_distance(x = centroid_points, y = WhitesVeg_Lines2)

WhitesEco_IntersectType <- WhitesEco_Intersect2 %>% select(GRIDCODE, WhiteEco_area)

Grid_10km2 = st_as_sf(Grid_10km)
SlaveExports_Nunn2 <- st_as_sf(SlaveExports_Nunn)
EWI_Intersect2 <- st_as_sf(EWI_Intersect)

Grid_10km2 <- st_join(Grid_10km2, centroid_points, join = st_intersects)
Grid_10km2 <- st_join(Grid_10km2, WhitesEco_IntersectType, join = st_intersects)
Grid_10km2 <- st_join(Grid_10km2, SlaveExports_Nunn2, join = st_intersects)
Grid_10km2 <- st_join(Grid_10km2, EWI_Intersect2, join = st_intersects)

save(Grid_10km2, file = "Data/Calculated Output/Grid_Covariates.RData") 

d <- as.data.frame(Grid_10km2) 
d <- d[ -c(78) ]
write.csv(d, "Data/Calculated Output/Covariates.csv", row.names = FALSE, na = "")

}

#Atlas
{
load("Data/Calculated Output/Grid_Covariates.RData")
  
  Atlas_6hr_Intersect2 = st_as_sf(Atlas_6hr_Intersect)
  Atlas_8hr_Intersect2 = st_as_sf(Atlas_8hr_Intersect)
  Atlas_10hr_Intersect2 = st_as_sf(Atlas_10hr_Intersect)
  Grid_10km2 = st_as_sf(Grid_10km)
  
  Atlas_6hr_Intersect2 <- Atlas_6hr_Intersect2 %>% select(Polity_6hr, Atlas_6hr_area_km, Atlas_6hr_area_km_Intersect)
  Atlas_8hr_Intersect2 <- Atlas_8hr_Intersect2 %>% select(Polity_8hr, Atlas_8hr_area_km, Atlas_8hr_area_km_Intersect)
  Atlas_10hr_Intersect2 <- Atlas_10hr_Intersect2 %>% select(Polity_10hr, Atlas_10hr_area_km, Atlas_10hr_area_km_Intersect)
  
  Grid_10km2 <- st_join(Grid_10km2, Atlas_6hr_Intersect2, join = st_intersects)
  Grid_10km2 <- st_join(Grid_10km2, Atlas_8hr_Intersect2, join = st_intersects)
  Grid_10km2 <- st_join(Grid_10km2, Atlas_10hr_Intersect2, join = st_intersects)
  
d <- as.data.frame(Grid_10km2) 
d <- d[ -c(16) ]
write.csv(d, "Data/Calculated Output/Atlas_Main.csv", row.names = FALSE, na = "")
}

#Atlas Composite
{
load("Data/Calculated Output/Grid_Covariates.RData")
  
sf_use_s2(FALSE)

Atlas_Comp_6hr_Intersect2 = st_as_sf(Atlas_Comp_6hr_Intersect)
Atlas_Comp_8hr_Intersect2 = st_as_sf(Atlas_Comp_8hr_Intersect)
Atlas_Comp_10hr_Intersect2 = st_as_sf(Atlas_Comp_10hr_Intersect)

Atlas_Comp_6hr_Intersect2 <- Atlas_Comp_6hr_Intersect2 %>% select(Atlas_Comp_6hr, Atlas_Comp_6hr_area_km, Atlas_Comp_6hr_area_km_Intersect)
Atlas_Comp_8hr_Intersect2 <- Atlas_Comp_8hr_Intersect2 %>% select(Atlas_Comp_8hr, Atlas_Comp_8hr_area_km, Atlas_Comp_8hr_area_km_Intersect)
Atlas_Comp_10hr_Intersect2 <- Atlas_Comp_10hr_Intersect2 %>% select(Atlas_Comp_10hr, Atlas_Comp_10hr_area_km, Atlas_Comp_10hr_area_km_Intersect)

Grid_10km2 <- st_join(Grid_10km2, Atlas_Comp_6hr_Intersect2, join = st_intersects)
Grid_10km2 <- st_join(Grid_10km2, Atlas_Comp_8hr_Intersect2, join = st_intersects)
Grid_10km2 <- st_join(Grid_10km2, Atlas_Comp_10hr_Intersect2, join = st_intersects)

d <- as.data.frame(Grid_10km2) 
d <- d[ -c(16) ]
write.csv(d, "Data/Calculated Output/Atlas_Composite.csv", row.names = FALSE, na = "")
}

#Atlas 1880
{
load("Data/Calculated Output/Grid_Covariates.RData")
  
sf_use_s2(FALSE)

Atlas_1880_6hr_Intersect2 = st_as_sf(Atlas_1880_6hr_Intersect)
Atlas_1880_8hr_Intersect2 = st_as_sf(Atlas_1880_8hr_Intersect)
Atlas_1880_10hr_Intersect2 = st_as_sf(Atlas_1880_10hr_Intersect)

Atlas_1880_6hr_Intersect2 <- Atlas_1880_6hr_Intersect2 %>% select(Atlas_1880_6hr, Atlas_1880_6hr_area_km, Atlas_1880_6hr_area_km_Intersect)
Atlas_1880_8hr_Intersect2 <- Atlas_1880_8hr_Intersect2 %>% select(Atlas_1880_8hr, Atlas_1880_8hr_area_km, Atlas_1880_8hr_area_km_Intersect)
Atlas_1880_10hr_Intersect2 <- Atlas_1880_10hr_Intersect2 %>% select(Atlas_1880_10hr, Atlas_1880_10hr_area_km, Atlas_1880_10hr_area_km_Intersect)

Grid_10km2 <- st_join(Grid_10km2, Atlas_1880_6hr_Intersect2, join = st_intersects)
Grid_10km2 <- st_join(Grid_10km2, Atlas_1880_8hr_Intersect2, join = st_intersects)
Grid_10km2 <- st_join(Grid_10km2, Atlas_1880_10hr_Intersect2, join = st_intersects)
Atlas_1880_6hr_Intersect2 <- Atlas_1880_6hr_Intersect2 %>% select(Atlas_1880_6hr, Atlas_1880_6hr_area_km, Atlas_1880_6hr_area_km_Intersect)
Atlas_1880_8hr_Intersect2 <- Atlas_1880_8hr_Intersect2 %>% select(Atlas_1880_8hr, Atlas_1880_8hr_area_km, Atlas_1880_8hr_area_km_Intersect)
Atlas_1880_10hr_Intersect2 <- Atlas_1880_10hr_Intersect2 %>% select(Atlas_1880_10hr, Atlas_1880_10hr_area_km, Atlas_1880_10hr_area_km_Intersect)

Grid_10km2 <- st_join(Grid_10km2, Atlas_1880_6hr_Intersect2, join = st_intersects)
Grid_10km2 <- st_join(Grid_10km2, Atlas_1880_8hr_Intersect2, join = st_intersects)
Grid_10km2 <- st_join(Grid_10km2, Atlas_1880_10hr_Intersect2, join = st_intersects)

d <- as.data.frame(Grid_10km2) 
d <- d[ -c(16) ]
write.csv(d, "Data/Calculated Output/Atlas_1880.csv", row.names = FALSE, na = "")
}

#Atlas Composite 1880
{
load("Data/Calculated Output/Grid_Covariates.RData")
  
sf_use_s2(FALSE)

Atlas_Comp_1880_6hr_Intersect2 = st_as_sf(Atlas_Comp_1880_6hr_Intersect)
Atlas_Comp_1880_8hr_Intersect2 = st_as_sf(Atlas_Comp_1880_8hr_Intersect)
Atlas_Comp_1880_10hr_Intersect2 = st_as_sf(Atlas_Comp_1880_10hr_Intersect)
Grid_10km2 = st_as_sf(Grid_10km)

Atlas_Comp_1880_6hr_Intersect2 <- Atlas_Comp_1880_6hr_Intersect2 %>% select(Atlas_Comp_1880_6hr, Atlas_Comp_1880_6hr_area_km, Atlas_Comp_1880_6hr_area_km_Intersect)
Atlas_Comp_1880_8hr_Intersect2 <- Atlas_Comp_1880_8hr_Intersect2 %>% select(Atlas_Comp_1880_8hr, Atlas_Comp_1880_8hr_area_km, Atlas_Comp_1880_8hr_area_km_Intersect)
Atlas_Comp_1880_10hr_Intersect2 <- Atlas_Comp_1880_10hr_Intersect2 %>% select(Atlas_Comp_1880_10hr, Atlas_Comp_1880_10hr_area_km, Atlas_Comp_1880_10hr_area_km_Intersect)

Grid_10km2 <- st_join(Grid_10km2, Atlas_Comp_1880_6hr_Intersect2, join = st_intersects)
Grid_10km2 <- st_join(Grid_10km2, Atlas_Comp_1880_8hr_Intersect2, join = st_intersects)
Grid_10km2 <- st_join(Grid_10km2, Atlas_Comp_1880_10hr_Intersect2, join = st_intersects)

d <- as.data.frame(Grid_10km2) 
d <- d[ -c(16) ]
write.csv(d, "Data/Calculated Output/Atlas_1880_Composite.csv", row.names = FALSE, na = "")
}

# Contested Territory
{
load("Data/Calculated Output/Grid_Covariates.RData")
  sf_use_s2(FALSE)
  
  Contested_1880_8hr_Intersect2 = st_as_sf(Contested_1880_8hr_Intersect)
  Contested_Composite_1880_8hr_Intersect2 = st_as_sf(Contested_Composite_1880_8hr_Intersect)

  Contested_1880_8hr_Intersect2 <- Contested_1880_8hr_Intersect2 %>% select(Polity_Contested_1880_8hr, Contested_1880_8hr_area_km_Intersect)
  Contested_Composite_1880_8hr_Intersect2 <- Contested_Composite_1880_8hr_Intersect2 %>% select(Polity_Contested_Comp_1880_8hr, Contested_Composite_1880_8hr_area_km_Intersect)

  Grid_10km2 <- st_join(Grid_10km2, Contested_1880_8hr_Intersect2, join = st_intersects)
  Grid_10km2 <- st_join(Grid_10km2, Contested_Composite_1880_8hr_Intersect2, join = st_intersects)

  d <- as.data.frame(Grid_10km2) 
  d <- d[ -c(44) ]
  write.csv(d, "Data/Calculated Output/Contested_Territory.csv", row.names = FALSE, na = "")
}

# Uncertainty Estimates
{
load("Data/Calculated Output/Grid_Covariates.RData")

  sf_use_s2(FALSE)
  
  UncertaintyEst_Max_Intersect2 = st_as_sf(UncertaintyEst_Max_Intersect)
  UncertaintyEst_Min_Intersect2 = st_as_sf(UncertaintyEst_Min_Intersect)
  
  UncertaintyEst_Max_Intersect2 <- UncertaintyEst_Max_Intersect2 %>% select(Polity_UncertaintyMax, UncertaintyEst_Max_area_km_Intersect)
  UncertaintyEst_Min_Intersect2 <- UncertaintyEst_Min_Intersect2 %>% select(Polity_UncertaintyMin, UncertaintyEst_Min_area_km_Intersect)
  
  Grid_10km2 <- st_join(Grid_10km2, UncertaintyEst_Max_Intersect2, join = st_intersects)
  Grid_10km2 <- st_join(Grid_10km2, UncertaintyEst_Min_Intersect2, join = st_intersects)
  
  d <- as.data.frame(Grid_10km2) 
  d <- d[ -c(44) ]
  write.csv(d, "Data/Calculated Output/Uncertainty_Estimates.csv", row.names = FALSE, na = "")
}
}

### IMPORT VIOLENCE DATA & MERGE TO GRID
#ACLED data
{
  load("Data/Calculated Output/Grid_Covariates.RData")
  
  ACLED_Base <- 
  read.csv("Data/Outcome Variables/Conflict Data/1997-01-01-2023-03-01-East_Asia-Eastern_Africa-Middle_Africa-Southern_Africa-Western_Africa_MW.csv",
           stringsAsFactors = FALSE)

  sf_use_s2(FALSE)
  
ACLED_Base_Plotted <- SpatialPointsDataFrame(ACLED_Base[,1:2], ACLED_Base)
crs(ACLED_Base_Plotted) <- CRS('+init=EPSG:4326')                                   
plot(ACLED_Base_Plotted)

ACLED_Base_Plotted2 = st_as_sf(ACLED_Base_Plotted)

ACLED_Base_Plotted2 %>% select(33) -> ACLED_Base_Plotted2

Grid_10km2_1 <- st_join(Grid_10km2, ACLED_Base_Plotted2, join = st_intersects)

d <- as.data.frame(Grid_10km2_1) 
d <- d[ -c(41) ]
write.csv(d, "Data/Calculated Output/AcledID.csv", row.names = FALSE, na = "")
}

#UCDP GED data
{
  load("Data/Calculated Output/Grid_Covariates.RData")
  
  GED_Base <- 
  read.csv("Data/Outcome Variables/Conflict Data/GED_23.1_Africa_Final_MW.csv",
           stringsAsFactors = FALSE)
GED_Base_Plotted <- SpatialPointsDataFrame(GED_Base[,1:2], GED_Base)
crs(GED_Base_Plotted) <- CRS('+init=EPSG:4326')                                   
plot(GED_Base_Plotted)

sf_use_s2(FALSE)

GED_Base_Plotted2 = st_as_sf(GED_Base_Plotted)

GED_Base_Plotted2 %>% select(3) -> GED_Base_Plotted2

Grid_10km2_2 <- st_join(Grid_10km2, GED_Base_Plotted2, join = st_intersects)

d <- as.data.frame(Grid_10km2_2) 
d <- d[ -c(74) ]
write.csv(d, "Data/Calculated Output/GEDID.csv", row.names = FALSE, na = "")
}
  
#META data
{
  load("Data/Calculated Output/Grid_Covariates.RData")
  
  Meta_Base <- 
      read.csv("Data/Outcome Variables/Meta Relative Wealth Data/Africa_Master_relative_wealth_index.csv",
               stringsAsFactors = FALSE)
    Meta_Base_Plotted <- SpatialPointsDataFrame(Meta_Base[,3:2], Meta_Base)
    crs(Meta_Base_Plotted) <- CRS('+init=EPSG:4326')                                   
    plot(Meta_Base_Plotted)
    
    sf_use_s2(FALSE)
    
    Meta_Base_Plotted2 = st_as_sf(Meta_Base_Plotted)
    
    Grid_10km2_2 <- st_join(Grid_10km2, Meta_Base_Plotted2, join = st_intersects)
    
    d <- as.data.frame(Grid_10km2_2) 
    d <- d[ -c(45) ]
    write.csv(d, "Data/Calculated Output/Meta_Wealth.csv", row.names = FALSE, na = "")
}

#AB data
{
  load("Data/Calculated Output/Grid_Covariates.RData")
  
  AB_clusters <- 
    read.csv("Data/Mechanisms/AB_Coordinates.csv",
             stringsAsFactors = FALSE)
  AB_clusters_Plotted <- SpatialPointsDataFrame(AB_clusters[,3:2], AB_clusters)
  crs(AB_clusters_Plotted) <- CRS('+init=EPSG:4326')                                   
  plot(AB_clusters_Plotted)
  
  sf_use_s2(FALSE)
  
  AB_clusters_Plotted2 = st_as_sf(AB_clusters_Plotted)
  Grid_10km2 = st_as_sf(Grid_10km)
  
  Grid_10km2_2 <- st_join(Grid_10km2, AB_clusters_Plotted2, join = st_intersects)
  
  d <- as.data.frame(Grid_10km2_2) 
  d <- d[ -c(10) ]
  write.csv(d, "Data/Calculated Output/AB_Cluster_Match.csv", row.names = FALSE, na = "")
  
  AB_Heaping <- 
    read.csv("Data/Mechanisms/AB_heaping.csv",
             stringsAsFactors = FALSE)
  AB_Heaping_Plotted <- SpatialPointsDataFrame(AB_Heaping[,4:3], AB_Heaping)
  crs(AB_Heaping_Plotted) <- CRS('+init=EPSG:4326')                                   
  plot(AB_Heaping_Plotted)
  
  sf_use_s2(FALSE)
  
  AB_Heaping_Plotted2 = st_as_sf(AB_Heaping_Plotted)
  
  Grid_10km2_2 <- st_join(Grid_10km2, AB_Heaping_Plotted2, join = st_intersects)
  
  d <- as.data.frame(Grid_10km2_2) 
  d <- d[ -c(83) ]
  write.csv(d, "Data/Calculated Output/AB_Heaping_Match.csv", row.names = FALSE, na = "")
  
AB_Heaping_Colonial <- 
    read.csv("Data/Mechanisms/AB_heaping_Colonial.csv",
             stringsAsFactors = FALSE)
  AB_Heaping_Colonial_Plotted <- SpatialPointsDataFrame(AB_Heaping_Colonial[,4:3], AB_Heaping_Colonial)
  crs(AB_Heaping_Colonial_Plotted) <- CRS('+init=EPSG:4326')                                   
  plot(AB_Heaping_Colonial_Plotted)
  
  sf_use_s2(FALSE)
  
  AB_Heaping_Colonial_Plotted2 = st_as_sf(AB_Heaping_Colonial_Plotted)
  
  Grid_10km2_2 <- st_join(Grid_10km2, AB_Heaping_Colonial_Plotted2, join = st_intersects)
  
  d <- as.data.frame(Grid_10km2_2) 
  d <- d[ -c(83) ]
  write.csv(d, "Data/Calculated Output/AB_Heaping_Colonial_Match.csv", row.names = FALSE, na = "")
  
}

#ADM1 (Mueller-Crepon)
{
ADM1 <- read_sf("Data/Grid Covariates/State Reach/admin_units/africa_regions_panel.GeoJSON")
year <- 2000
regcap2000 <- raster("Data/Grid Covariates/State Reach/state_reach/time2regcap_1945_2016_dynroads.tif",
                     band =  year-1945+1)
natcap2000 <- raster("Data/Grid Covariates/State Reach/state_reach/time2natcap_1946_2016_dynroads.tif",
                     band =  year-1945+1)
year <- 1966
regcap1966 <- raster("Data/Grid Covariates/State Reach/state_reach/time2regcap_1945_2016_dynroads.tif",
                     band =  year-1945+1)
natcap1966 <- raster("Data/Grid Covariates/State Reach/state_reach/time2natcap_1946_2016_dynroads.tif",
                     band =  year-1945+1)
year <- 1977
regcap1977 <- raster("Data/Grid Covariates/State Reach/state_reach/time2regcap_1945_2016_dynroads.tif",
                     band =  year-1945+1)
natcap1977 <- raster("Data/Grid Covariates/State Reach/state_reach/time2natcap_1946_2016_dynroads.tif",
                     band =  year-1945+1)

year <- 1980
regcap1980 <- raster("Data/Grid Covariates/State Reach/state_reach/time2regcap_1945_2016_dynroads.tif",
                     band =  year-1945+1)
natcap1980 <- raster("Data/Grid Covariates/State Reach/state_reach/time2natcap_1946_2016_dynroads.tif",
                     band =  year-1945+1)
year <- 1968
regcap1968 <- raster("Data/Grid Covariates/State Reach/state_reach/time2regcap_1945_2016_dynroads.tif",
                     band =  year-1945+1)
natcap1968 <- raster("Data/Grid Covariates/State Reach/state_reach/time2natcap_1946_2016_dynroads.tif",
                     band =  year-1945+1)
year <- 1975
regcap1975 <- raster("Data/Grid Covariates/State Reach/state_reach/time2regcap_1945_2016_dynroads.tif",
                     band =  year-1945+1)
natcap1975 <- raster("Data/Grid Covariates/State Reach/state_reach/time2natcap_1946_2016_dynroads.tif",
                     band =  year-1945+1)

Grid_10km$regcap1975 <- exact_extract(regcap1975,Grid_10km,fun="mean")
Grid_10km$natcap1975 <- exact_extract(natcap1975,Grid_10km,fun="mean")
Grid_10km$regcap1968 <- exact_extract(regcap1968,Grid_10km,fun="mean")
Grid_10km$natcap1968 <- exact_extract(natcap1968,Grid_10km,fun="mean")
Grid_10km$regcap1980 <- exact_extract(regcap1980,Grid_10km,fun="mean")
Grid_10km$natcap1980 <- exact_extract(natcap1980,Grid_10km,fun="mean")
Grid_10km$regcap1977 <- exact_extract(regcap1977,Grid_10km,fun="mean")
Grid_10km$natcap1977 <- exact_extract(natcap1977,Grid_10km,fun="mean")
Grid_10km$regcap1966 <- exact_extract(regcap1966,Grid_10km,fun="mean")
Grid_10km$natcap1966 <- exact_extract(natcap1966,Grid_10km,fun="mean")

sf_use_s2(FALSE)
Grid_10km_admin = st_as_sf(Grid_10km)

d <- as.data.frame(Grid_10km_admin) 
d <- d[ -c(3, 14) ]
write.csv(d, "Data/Calculated Output/Grid_StateReach.csv", row.names = FALSE, na = "")


ADM1 <- shapefile("Data/Grid Covariates/State Reach/admin_units/Admin_2016_alt.shp")
 
sf_use_s2(FALSE)
ADM1_2 = st_as_sf(ADM1)

ADM1_centroid <- st_join( centroid_points, ADM1_2, join = st_intersects)

Grid_10km2 = st_as_sf(Grid_10km)

Grid_10km2_2 <- st_join(Grid_10km2, ADM1_centroid, join = st_intersects)

d <- as.data.frame(Grid_10km2_2) 
d <- d[ -c(15) ]
write.csv(d, "Data/Calculated Output/Grid_Admin_Match_ACPED.csv", row.names = FALSE, na = "")

}

#Horse data
{
  load("Data/Calculated Output/Grid_Covariates.RData")
  
  sf_use_s2(FALSE)
  
  Blecher_Pony <- shapefile("Data/Grid Covariates/Horses/Blecher_1993_Horses.shp")
  Pony_Intersect <- intersect(Blecher_Pony, Grid_10km)
  Pony_Intersect2 = st_as_sf(Pony_Intersect)
  Pony_Intersect2 <- Pony_Intersect2 %>% select(Horses)

  Grid_10km2 = st_as_sf(Grid_10km)

  Grid_10km2 <- st_join(Grid_10km2, Pony_Intersect2, join = st_intersects)
  
  d <- as.data.frame(Grid_10km2) 
  d <- d[ -c(6) ]
  write.csv(d, "Data/Calculated Output/Horses_Grid_Overlap.csv", row.names = FALSE, na = "")
}

